satkit 0.16.2

Satellite Toolkit
Documentation
{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "# Lambert Targeting"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1",
   "metadata": {},
   "source": [
    "This tutorial demonstrates how to use satkit's Lambert solver to design orbit transfers.\n",
    "We'll compute the delta-v required for a satellite to move between orbits and visualize the transfer geometry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import scienceplots  # noqa: F401\n",
    "plt.style.use([\"science\", \"no-latex\", \"../satkit.mplstyle\"])\n",
    "%config InlineBackend.figure_formats = ['svg']\n",
    "\n",
    "import satkit as sk\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3",
   "metadata": {},
   "source": [
    "## Scenario: LEO to MEO Transfer\n",
    "\n",
    "A satellite is in a 400 km circular orbit (LEO) and needs to transfer to a 2000 km circular orbit (MEO).\n",
    "We'll compute the optimal transfer for a range of transfer angles and find the one that minimizes total delta-v."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Orbit parameters\n",
    "alt_1 = 400e3    # departure altitude (m)\n",
    "alt_2 = 2000e3   # arrival altitude (m)\n",
    "r_earth = sk.consts.earth_radius\n",
    "mu = sk.consts.mu_earth\n",
    "\n",
    "r1_mag = r_earth + alt_1\n",
    "r2_mag = r_earth + alt_2\n",
    "\n",
    "# Circular velocities at each orbit\n",
    "v_circ_1 = np.sqrt(mu / r1_mag)\n",
    "v_circ_2 = np.sqrt(mu / r2_mag)\n",
    "\n",
    "print(f\"Departure orbit: {alt_1/1e3:.0f} km altitude, v = {v_circ_1:.1f} m/s\")\n",
    "print(f\"Arrival orbit:   {alt_2/1e3:.0f} km altitude, v = {v_circ_2:.1f} m/s\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5",
   "metadata": {},
   "source": [
    "## Delta-v vs Transfer Angle\n",
    "\n",
    "We place the satellite at a fixed departure point and sweep the arrival point around the higher orbit,\n",
    "computing the Lambert solution and delta-v for each transfer angle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Departure position (on x-axis)\n",
    "r1 = np.array([r1_mag, 0.0, 0.0])\n",
    "\n",
    "# Sweep transfer angles from 30 to 170 degrees\n",
    "angles_deg = np.linspace(30, 180, 200)\n",
    "angles_rad = np.radians(angles_deg)\n",
    "\n",
    "# Hohmann transfer time (reference)\n",
    "a_transfer = (r1_mag + r2_mag) / 2\n",
    "t_hohmann = np.pi * np.sqrt(a_transfer**3 / mu)\n",
    "\n",
    "dv1_arr = []\n",
    "dv2_arr = []\n",
    "dv_total_arr = []\n",
    "\n",
    "for theta in angles_rad:\n",
    "    r2 = np.array([r2_mag * np.cos(theta), r2_mag * np.sin(theta), 0.0])\n",
    "\n",
    "    # Scale TOF with transfer angle: fraction of Hohmann time\n",
    "    tof = t_hohmann * (theta / np.pi)\n",
    "\n",
    "    solutions = sk.lambert(r1, r2, tof)\n",
    "    v1_transfer, v2_transfer = solutions[0]\n",
    "\n",
    "    # Departure: satellite in circular orbit, moving in +y direction\n",
    "    v1_circular = np.array([0.0, v_circ_1, 0.0])\n",
    "\n",
    "    # Arrival: circular velocity tangent to orbit at arrival point\n",
    "    v2_circular = np.array([\n",
    "        -v_circ_2 * np.sin(theta),\n",
    "        v_circ_2 * np.cos(theta),\n",
    "        0.0,\n",
    "    ])\n",
    "\n",
    "    dv1 = np.linalg.norm(v1_transfer - v1_circular)\n",
    "    dv2 = np.linalg.norm(v2_transfer - v2_circular)\n",
    "\n",
    "    dv1_arr.append(dv1)\n",
    "    dv2_arr.append(dv2)\n",
    "    dv_total_arr.append(dv1 + dv2)\n",
    "\n",
    "dv1_arr = np.array(dv1_arr)\n",
    "dv2_arr = np.array(dv2_arr)\n",
    "dv_total_arr = np.array(dv_total_arr)\n",
    "\n",
    "# Find minimum total delta-v\n",
    "idx_min = np.argmin(dv_total_arr)\n",
    "print(f\"Minimum total delta-v: {dv_total_arr[idx_min]:.1f} m/s at {angles_deg[idx_min]:.1f} degrees\")\n",
    "print(f\"  Departure burn: {dv1_arr[idx_min]:.1f} m/s\")\n",
    "print(f\"  Arrival burn:   {dv2_arr[idx_min]:.1f} m/s\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 5))\n",
    "ax.plot(angles_deg, dv1_arr, label=r\"$\\Delta v_1$ (departure)\")\n",
    "ax.plot(angles_deg, dv2_arr, label=r\"$\\Delta v_2$ (arrival)\")\n",
    "ax.plot(angles_deg, dv_total_arr, label=r\"$\\Delta v_{total}$\", linewidth=2)\n",
    "ax.axvline(angles_deg[idx_min], color=\"#BBBBBB\", linestyle=\":\", linewidth=1)\n",
    "ax.set_xlabel(\"Transfer Angle (degrees)\")\n",
    "ax.set_ylabel(r\"$\\Delta v$ (m/s)\")\n",
    "ax.set_title(\"LEO to MEO Transfer\")\n",
    "ax.legend()\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8",
   "metadata": {},
   "source": [
    "## Pork-Chop Plot: Delta-v vs Transfer Angle and Time\n",
    "\n",
    "A **pork-chop plot** shows total delta-v as a function of both transfer angle and time of flight,\n",
    "revealing the trade-off between fuel cost and transfer duration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sweep both transfer angle and time of flight\n",
    "angles = np.linspace(60, 300, 120)\n",
    "tofs = np.linspace(0.4 * t_hohmann, 3.0 * t_hohmann, 100)\n",
    "\n",
    "DV = np.full((len(tofs), len(angles)), np.nan)\n",
    "\n",
    "for i, t in enumerate(tofs):\n",
    "    for j, theta_deg in enumerate(angles):\n",
    "        theta = np.radians(theta_deg)\n",
    "        r2 = np.array([r2_mag * np.cos(theta), r2_mag * np.sin(theta), 0.0])\n",
    "\n",
    "        try:\n",
    "            solutions = sk.lambert(r1, r2, t)\n",
    "            v1_t, v2_t = solutions[0]\n",
    "\n",
    "            v1_c = np.array([0.0, v_circ_1, 0.0])\n",
    "            v2_c = np.array([\n",
    "                -v_circ_2 * np.sin(theta),\n",
    "                v_circ_2 * np.cos(theta),\n",
    "                0.0,\n",
    "            ])\n",
    "\n",
    "            dv = np.linalg.norm(v1_t - v1_c) + np.linalg.norm(v2_t - v2_c)\n",
    "            if dv < 15000:\n",
    "                DV[i, j] = dv\n",
    "        except Exception:\n",
    "            pass\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(9, 6))\n",
    "levels = np.linspace(500, 8000, 20)\n",
    "cs = ax.contourf(angles, tofs / 60, DV, levels=levels, cmap=\"viridis\", extend=\"max\")\n",
    "ax.contour(angles, tofs / 60, DV, levels=levels, colors=\"white\", linewidths=0.3)\n",
    "cbar = fig.colorbar(cs, ax=ax, label=r\"Total $\\Delta v$ (m/s)\")\n",
    "ax.set_xlabel(\"Transfer Angle (degrees)\")\n",
    "ax.set_ylabel(\"Time of Flight (minutes)\")\n",
    "ax.set_title(\"LEO to MEO Transfer: Pork-Chop Plot\")\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10",
   "metadata": {},
   "source": [
    "## Transfer Orbit Visualization\n",
    "\n",
    "Let's visualize the optimal transfer orbit alongside the departure and arrival orbits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Optimal transfer\n",
    "theta_opt = np.radians(angles_deg[idx_min])\n",
    "r2_opt = np.array([r2_mag * np.cos(theta_opt), r2_mag * np.sin(theta_opt), 0.0])\n",
    "\n",
    "tof_opt = t_hohmann * (theta_opt / np.pi)\n",
    "solutions = sk.lambert(r1, r2_opt, tof_opt)\n",
    "v1_transfer, v2_transfer = solutions[0]\n",
    "\n",
    "# Propagate transfer orbit via two-body integration\n",
    "from scipy.integrate import solve_ivp\n",
    "\n",
    "def twobody(t, state):\n",
    "    r = state[:3]\n",
    "    v = state[3:]\n",
    "    r_norm = np.linalg.norm(r)\n",
    "    a = -mu / r_norm**3 * r\n",
    "    return np.concatenate([v, a])\n",
    "\n",
    "state0 = np.concatenate([r1, v1_transfer])\n",
    "sol = solve_ivp(twobody, [0, tof_opt], state0, rtol=1e-10, atol=1e-12,\n",
    "                dense_output=True)\n",
    "t_plot = np.linspace(0, tof_opt, 500)\n",
    "positions = sol.sol(t_plot)[:3].T\n",
    "\n",
    "# Draw orbits\n",
    "fig, ax = plt.subplots(figsize=(8, 8))\n",
    "theta_circ = np.linspace(0, 2 * np.pi, 500)\n",
    "\n",
    "# Earth\n",
    "earth = plt.Circle((0, 0), r_earth / 1e3, color=\"lightblue\", ec=\"steelblue\", linewidth=1)\n",
    "ax.add_patch(earth)\n",
    "\n",
    "# Departure orbit\n",
    "ax.plot(r1_mag * np.cos(theta_circ) / 1e3, r1_mag * np.sin(theta_circ) / 1e3,\n",
    "        \"--\", color=\"#BBBBBB\", linewidth=1, label=f\"LEO ({alt_1/1e3:.0f} km)\")\n",
    "\n",
    "# Arrival orbit\n",
    "ax.plot(r2_mag * np.cos(theta_circ) / 1e3, r2_mag * np.sin(theta_circ) / 1e3,\n",
    "        \"--\", color=\"#BBBBBB\", linewidth=1, label=f\"MEO ({alt_2/1e3:.0f} km)\")\n",
    "\n",
    "# Transfer arc\n",
    "ax.plot(positions[:, 0] / 1e3, positions[:, 1] / 1e3,\n",
    "        linewidth=2, color=\"#CC3311\", label=\"Transfer orbit\")\n",
    "\n",
    "# Endpoints\n",
    "ax.plot(r1[0] / 1e3, r1[1] / 1e3, \"o\", color=\"#0077BB\", markersize=8, zorder=5)\n",
    "ax.annotate(\"Departure\", xy=(r1[0] / 1e3, r1[1] / 1e3),\n",
    "            xytext=(15, 15), textcoords=\"offset points\", fontsize=11)\n",
    "ax.plot(r2_opt[0] / 1e3, r2_opt[1] / 1e3, \"s\", color=\"#009988\", markersize=8, zorder=5)\n",
    "ax.annotate(\"Arrival\", xy=(r2_opt[0] / 1e3, r2_opt[1] / 1e3),\n",
    "            xytext=(15, -20), textcoords=\"offset points\", fontsize=11)\n",
    "\n",
    "ax.set_xlabel(\"x (km)\")\n",
    "ax.set_ylabel(\"y (km)\")\n",
    "ax.set_title(f\"Transfer Orbit  ({angles_deg[idx_min]:.0f}°, \" + r\"$\\Delta v$ = \" + f\"{dv_total_arr[idx_min]:.0f} m/s)\")\n",
    "ax.set_aspect(\"equal\")\n",
    "ax.legend(loc=\"lower left\")\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12",
   "metadata": {},
   "source": [
    "## Interplanetary Example: Earth to Mars\n",
    "\n",
    "Lambert targeting also works for interplanetary transfers by using the Sun's gravitational parameter and heliocentric positions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13",
   "metadata": {},
   "outputs": [],
   "source": [
    "mu_sun = sk.consts.mu_sun\n",
    "au = sk.consts.au\n",
    "\n",
    "# Simplified circular coplanar orbits\n",
    "r_earth_orbit = 1.0 * au\n",
    "r_mars_orbit = 1.524 * au\n",
    "\n",
    "# Earth at departure\n",
    "r1_helio = np.array([r_earth_orbit, 0, 0])\n",
    "\n",
    "# Mars at arrival (~135 degree transfer)\n",
    "transfer_angle = np.radians(135)\n",
    "r2_helio = np.array([\n",
    "    r_mars_orbit * np.cos(transfer_angle),\n",
    "    r_mars_orbit * np.sin(transfer_angle),\n",
    "    0,\n",
    "])\n",
    "\n",
    "# Transfer time: ~8 months\n",
    "tof_days = 245\n",
    "tof_interp = tof_days * 86400\n",
    "\n",
    "solutions = sk.lambert(r1_helio, r2_helio, tof_interp, mu=mu_sun)\n",
    "v1_t, v2_t = solutions[0]\n",
    "\n",
    "# Circular velocities\n",
    "v_earth = np.sqrt(mu_sun / r_earth_orbit)\n",
    "v_mars = np.sqrt(mu_sun / r_mars_orbit)\n",
    "\n",
    "v1_earth = np.array([0, v_earth, 0])\n",
    "v2_mars = np.array([\n",
    "    -v_mars * np.sin(transfer_angle),\n",
    "    v_mars * np.cos(transfer_angle),\n",
    "    0,\n",
    "])\n",
    "\n",
    "dv_depart = np.linalg.norm(v1_t - v1_earth)\n",
    "dv_arrive = np.linalg.norm(v2_t - v2_mars)\n",
    "\n",
    "print(f\"Earth-Mars transfer ({tof_days} days):\")\n",
    "print(f\"  Departure delta-v:  {dv_depart/1e3:.2f} km/s\")\n",
    "print(f\"  Arrival delta-v:    {dv_arrive/1e3:.2f} km/s\")\n",
    "print(f\"  Total delta-v:      {(dv_depart + dv_arrive)/1e3:.2f} km/s\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.12.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}